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ABSTRACT 

This study concerns the determination of the contact stresses and contact region around bolt holes and the 
bolt load distribution in single- and double-lap joints of composite laminates with arbitrarily located bolts 
under general mechanical loading conditions and uniform temperature change. The unknown contact 
stress distribution and contact region between the bolt and laminates and the interaction among the bolts 
require the bolt load distribution, as well as the contact stresses, to be as part of the solution. The present 
method is based on the complex potential theory and the variational formulation in order to account for 
bolt stiffness, bolt-hole clearance, and finite geometry of the composite laminates. 

1. INTRODUCTION 

Bolted joints provide the primary means for transferring load among composite components in the 
construction of aircraft and space structures. The stress state in a bolted joint is dependent primarily on 
the dimensions of the planar geometry, loading conditions, degree of material anisotropy, bolt-hole 
clearance, bolt flexibility, and friction between the laminates. Also, aircraft and space vehicles traveling at 
supersonic and hypersonic speeds can experience high temperature excursions. The influence of thermal 
expansions can be significant and may differ significantly among the materials for the bolts and 
laminates. As a result, high thermal stresses may develop as the temperature increases and may alter the 
bolt load distribution. Therefore, accurate determination of the stresses in bolted laminates under both 
mechanical and thermal loading is essential for reliable strength evaluation and failure prediction. 

A considerable amount of work on the behavior of composite joints with a single bolt exists in the 
literature. These studies investigated the stress distribution around a pin-loaded hole in laminated 
composites based on either finite element analysis or analytical methods. Since the contact stress 
distribution and the contact region are not known a priori , a majority of the models did not directly 
impose the boundary conditions appropriate for modeling the contact and non-contact regions between 
the bolt and the boundary of the hole. These models usually assumed a cosinusoidal bearing stress 
distribution or zero radial displacements over the contact region of the hole boundary. In the case of 
multi-bolt joints, the commonly accepted approach is to first determine the load distribution among the 
bolts in order to identify the critical (most highly loaded) bolt for a subsequent single-bolt analysis for 
local stress distribution. However, this type of analysis disregards the interaction among the bolts located 
in close proximity to each other. In order to eliminate these shortcomings, Madenci et al. (1998) 
developed a method for single-lap joints based on the boundary collocation technique. Their method 
determines the contact stresses and contact region, as well as the bolt load distribution, as part of the 
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solution procedure. However, this method fails to provide converged solutions consistently depending on 
the number of bolts and their location in relation to each other or to the free boundaries. A detailed 
validation and demonstration of their approach, as well as an extensive review of previous analyses, were 
reported in detail by Madenci et al. (1997). 

In the literature, there are essentially no direct analyses of double-lap bolted joints for solid laminates 
under general loading conditions and appropriate boundary conditions arising from contact phenomenon. 
Madenci et al. (1999) extended their boundary collocation technique for single-lap joints to consider 
double-lap joints and thermal loading. This method provided converged results for particular 
configurations, but also suffered from consistent convergence arising from the explicit partitioning of the 
domain. 

Xiong and Poon (1994) introduced an analytical approach utilizing a variational formulation in 
conjunction with the complex potential theory to single- and double-lap joints with many bolts. Their 
approach considers each laminate of the joint separately. The coupling of the laminates is achieved 
through bolt displacements, which are permitted only in the direction of loading. In their two-stage 
analysis, the first stage provides the local deformation along the hole boundaries of one of the laminates 
subjected to the external boundary conditions and the prescribed cosinusoidal bearing stress representing 
the bolt load at each hole boundary. The local deformations and the bolt deflections are imposed as 
displacement constraints in the subsequent second stage to determine the contact stresses (bolt loads) and 
the contact region in the second laminate. Subsequently, these fastener loads are imposed as prescribed 
cosinusoidal bearing stress for the first stage of the analysis, and the iterative process continues until the 
constraint conditions are satisfied. 

This study presents an analysis method for determining the bolt load distribution in single- and 
double-lap joints while accounting for the contact phenomenon and the interaction among the bolts 
explicitly under bearing and by-pass loading with or without thermal loading. It is an extension of the 
analysis introduced by Xiong and Poon (1994) and eliminates the requirement of a two-stage analysis 
and the associated iterative process. The resulting equations are solved in a coupled manner, leading to 
the contact stresses, contact region, and bolt load distribution. 

2. PROBLEM STATEMENT 

The geometry of a bolted single- and double-lap joint with composite laminates is described in Fig. 1. The 
joint can be subjected to a combination of bearing, by-pass and shear loads, and a uniform temperature 
change. Each laminate of the single- and double-lap joints, joined with L number of bolts, is subjected to 
traction components, ( a = x,y ) along the segment of the external boundary of each region 
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denoted by r^. The section of the external boundary subjected to displacement constraints, it f 
( p = n,s ), is denoted by f (/:) . The subscripts “n" and “s” denote directions outward normal and 
tangent to the boundary, respectively. Each region with an area of can be under uniform 

temperature change, . The thickness of the laminates (regions) is denoted by hr ' . The contact 
region between the f th bolt and the hole boundary in the A: 111 region is denoted by . The sub- or 
superscripts “(it)” and “(()” refer to the regions (laminates) and bolts, respectively. Their ranges are 
specified by it = 1,..., AT and t = 1,..., L , with K and L being the total number of regions (laminates) and 
bolts, respectively. As illustrated in Fig. 2, the hole radius, a (, which is slightly larger than the bolt 
radius, R ( , leads to a clearance of S f . The hole and bolt radii remain the same for each region. As shown 
in Fig. 2, the center of each hole, located at (X f ,F f ), coincides with the origin of the Cartesian 
coordinates (jty,y^). 




Fig. 2. Position of a bolt before and after the load is executed. 

As shown in Fig. 3, the free-body diagram of each component of a lap joint, the unknown boundary 
traction components, If , arise from the deformation of the boundary given by - (u^ -w^) 
along a portion of the external boundary, f ^ . The unknown traction component in the outward normal 
direction, , arises from the deformation of the contact zone between the £ th bolt and the hole 
boundary in the it th region laminate. This contact zone deformation is expressed by = 
u n k>) along the contact region denoted by f ^ . The extent of the contact region 

denoted by is dependent on the bolt displacement, , deformation of the hole boundary, 
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Fig. 3. Free-body diagram of each component in a bolted single-lap joint. 

and the gap, S*j. () . Because of the absence of friction between the bolt and the laminate, the 
tangential component of the bolt displacement, , and the traction vector, , vanish, i.e. = 0 
and Xj*^=0. As shown in Fig. 4, at the point of initial contact (prior to any deformation of hole 
boundary), the gap between the hole boundary and the bolt (distance PP') in the k 0 ' region is defined by 

~PP' = S* k() =S ( [ l-cos(0-0*)] (1) 
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Fig. 4. The gap between the bolt and hole boundary immediately before and after the load is exerted 


in which 9* specifies the line of action, m , and 8^ is the clearance. The extent of the contact region, 
f , is defined by the angles d A and d B . The flexible bolts experience deflections given by 




x\ ’ *2 

A (0 r _. A (0 A (0 # A (0 A (0 A (0 ) 

-i^ x ] ’^*2’ *3’ vl > a v2’ a v3' 


yl 

y\ ’^ y 2'“y3 


for a single-lap joint 
for a double-lap joint 

,th 


(2a) 

(2b) 


with A^ } and A^ (i = 1, K) representing the bolt deflection components at the i ul point along the 
length of the t * bolt along the x- and y-directions, respectively. 

The material properties of each laminate are represented by the matrix relating the stress 

resultants, , to strain resultants, , with a,J3 = x,y , in the form 
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where Ajj are the components of the in-plane stiffness matrix A'*' of the ft* region. The strain 
components arising from temperature change, £' a p , are expressed as 
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where the coefficients a ( J ^ , a ^) , and <ar~ ; represent the thermal expansions of the ft m region with 
respect to the global (X ,Y) coordinate system. The corresponding thermal stress resultants are defined as 
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In matrix form, this relationship is expressed as 


V»=A ( « V*> 


( 6 ) 


in which *e^ = with } . 
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The stiffness matrix of the bolts, , is given by 


b^=l 


} (0 

J x 

0 


JO 


( 7 ) 


whose coefficients are determined by modeling the bolt as a beam under concentrated forces. The explicit 
expressions for bolt stiffness for a single- and double-lap joint, as well as the general lap configurations, 
are derived in Appendix A. 

These angles, the contact stresses, the components of bolt displacement, and the forces exerted by the 
bolts are the unknowns to be determined as part of the solution. Unless indicated otherwise, the subscripts 
a and j5 vary as a,fl = x,y, representing the (x,y) global coordinates. The subscript p varies as 
p = n,s , representing the directions normal and tangent to the boundary, as shown in Fig. 1. Also, only 
repeated subscripts imply summation. 


3. SOLUTION METHOD 

The solution method is based on the variational formulation in conjunction with the complex potential 
theory. The governing equations are derived by requiring the first variation of the total potential energy, 
arising from thermal and mechanical loads, to vanish. The in-plane equilibrium equations in each region 
are satisfied exactly by employing complex potential functions in the form suggested by Lekhnitskii 
(1968). However, each of the bolt equilibrium equations and the boundary conditions are satisfied by 
minimizing the total potential energy. 


3.1. Governing Equations 

The total potential energy for K regions connected with L number of bolts in the absence of friction 
between the laminates and between the laminates and the bolts along the contact region under mechanical 
and thermal loading can be expressed as 

K L K L K K 

k= 1 ^=1 k =\£= 1 k = I k = 1 

The strain energy of the k ^ laminate, , is given by 

U (k) ={ I N aP e aP dA ~ \ N a/S* £ afi dA < 9 > 




(*) 


( 8 ) 
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with its first variation as derived in Appendix B, 


su (k) =- J N$ j3 Su ( a k) dA+ J */<?>&!<?> dr 

A (k) f(*) 

L (k) 

+ J ^SufdT+Y, J 'tfSufdT 

p(A:) i ~ I p(^0 

in which a,fl = x,y and p = n,s . The strain energy of the l Ul bolt, , is given by 


with its first variation 


B<'>= 

2 ' ,J J 


with i,j = \,2K 


SB^=b\pAfS^ 


( 10 ) 


( 11 ) 


( 12 ) 


The potential of the reaction forces, AjfO and ( p = n,s ) , arising from the contact between the bolt 

and the hole boundary and the applied displacement constraints along the external boundary are denoted 
by and , respectively. They are expressed in the form 


in which 


and 




f(«> 



j(k() m 


contact between the fc th plate and ^ th bolt 
no contact between the & th plate and £ ih bolt 


= J dT with p = n,s 

f (*) 


(13) 


(14) 


The potential of the externally applied tractions, t ^ , denoted by is expressed as 

W^=- | t^u^dT with a-x,y 


(15) 


"(*) 


in which and u ^ represent the applied traction and displacement components in the x- and y- 
directions, respectively. Their first variations are obtained as 
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(16) 


<*>“«= J ;<«){»«) -i<‘o (4 W)_ (5 f Jlo ]4S*' ) < / r 
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and 


SW (k) — J t^ k) Su {k) dT (18) 

p(^) 

The first variation of the total potential energy can be obtained as 
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Noting that 5u \ k) , Su^, Su^„ k \ 3u[ k \ SA and SA ( n kf:) are arbitrary independent quantities and 
requiring the first variation of the total potential energy to vanish lead to the equilibrium equations for the 
laminates and bolts as 


N ik) -0 
N ccp,p ~ U 


on A 


(A) 


(20) 
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( 21 ) 


Final Report NAG-1 -2052 


Page 10 


Their associated boundary conditions are obtained as 

{* 4 * > -' 3‘ ) }=0 onP <‘> ( 22 a ) 

= Q on f (22b) 

*,<‘>=0 onf (W> ( 22 c ) 

{%* , +$ > }-0 onf <« ( 22 d ) 

{»«>-C«>}=0o n f<‘> (22c) 

-i?')(Af )+*('„,] =0 on f<«> (22f) 

where a, (3 = x>y\ p = n,s\ k=l,K; and £ = 1,...,L. 


3.2. Total Potential Energy 

The strain energy expression given in Eq. (9) can be rewritten in terms of displacement components as 


| *$■$*<“- f 


Its integration by parts yields 


t' , ‘ ) / [W? - 2 *^ ) K > 


<£4 


A<*> 


■/? 


Under uniform temperature distribution and applying Gauss’ theorem, it reduces to 

U W =~ J f ( N $ 


i(k) 


r(k) 


(23) 


(24) 


(25) 


This expression can be further simplified by invoking the stress resultants, N a a g, that satisfy the 
equilibrium equations, N a p p - 0 , given by Eq. (20), as 


U 


(*) 


4 1 

Y^k) 


(26) 
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The explicit expressions for the stress resultants and displacement components satisfying both the 
equilibrium equations and the compatibility conditions are given in Appendix C. 

Finally, the total potential energy expression given in Eq. (19) is reduced to a form, free of area 
integrals, as 


K 


N 


J r + l£if A<»Af +1 J 

*=i z r (*) L e=\ k=\f(k) 

*=1^=1 f(*£) *=If(*) 

with a,fi = x,y, p = n,s; k and i = 1,...,L . It can be rewritten in matrix form as 

L 
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where 
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and the matrix of unit normals, n 1 ^ , is expressed as 
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(29a) 

(29b) 

(29c) 

(29d) 


(29e) 


As given by Equation (3), the unknown bolt displacements are contained in vector and matrix b ( ^ 
represents the bolt stiffness. 
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The constraint condition given by Equation (22e), 


dp'* =Up -Up^ =0 with p = n,s 


(30) 


is rewritten in vector form as 


in which 


c <*>= u '(*>- u ' ( *>=0 


„'(*) =T (*) U (*) 


(31) 


(32) 


with T^- 1 representing the transformation matrix between the (x, y) and (n,s) coordinate systems. The 
displacement vectors u'^ and are defined in the (n,s) and (x,y) coordinate systems, respectively, 




(33) 


u<‘> r ] 


(34) 


Substituting for from Equation (B31) into Equation (32) yields 

u '(*) =T (*) u (*) 7 a (*) 


(35) 


leading to 


f (*) =x (*) u (*) r a (t) -u' w 


(36) 


The unknown traction vector, X (k> , defined in Eq. (29c) can be assumed as 

k ( *> = £ VjAf 

7=0 


(37) 


- /h\ 

where the matrix P j and vector of unknown coefficients A j are defined as 


P 7 = 


Pj 0 
0 Pj 


and \f )T ={A nj ,A sj } 


(38) 


with pj being the y^-order Legendre polynomial. In matrix notation, this equation becomes 
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where 


i (k) =PA (k) 


(39) 


P = [P 0 P, P 2 P j] and A ( * )T ={a£, A[, A£ Aj J (40) 

From Eqs. (36) and (40), the boundary integral of the product that appears in the expression for 

the total potential energy is obtained as 

Ji« )r c (i Vr= | A<‘> I 'p r T< i 'u ( ‘> 7 ’a«>jr- j A«> r p r T <l) a ( ‘Vr 

f(*) f (*) f (*) (41) 

= A (k)T C (k) a (k) -A (k)r f (k) 

where 

c (k) = | P r T (A:) U ( * >r </r (42) 

ft*) 

and 

f (k) = J P r T (k) u (k) dr ( 43 ) 

p(k) 


The constraint condition given by Equation (22f), 

=«?)-«£“>( ^P) +s m=o 

is rewritten in vector form as 


e («) =u ^)_ a ^)( A (O) + 5 ^ )=0 


(44) 


(45) 


The bolt displacement vector, , at the contact region can be expressed in terms of the bolt 

deflections in the form 


u' ( *°=G (W) A (0 


(46) 


in which the matrix G v and the vector A v J are defined as 
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G {k{) = 


sm 6> k [ 

cos 9 k i 

0 

0 


with J](Zj - z k ) defined as 


fj(zi~z k ) rj(z 2 -z k ) rj(z K -z k ) 

0 0 ••• 0 
0 0 0 

V(z\~z k ) r)(z 2 -z k ) ••• rj(z K -z k )_ 


f 1 if i = k 

0=1 K) 


and 

AO T -[AO AO ... AO AO AO ... AO] 

A _ \ A ;tl A *2 a xK a yl A >2 A >*:j 

With substitutions from Eq. (35) and (46), the expression for becomes 

S («) =T <*> v m T a m +«,*„, 

The unknown traction vector, , defined in Eq. (29c) can be assumed as 

l 


x«'>=X F i A ! 

i=0 


m 


* (bf\ 

where the matrix F,- and vector of unknown coefficients AJ ; are defined as 

and k\ k()T ={A^ £) ,o} 


F f = 


F t 0 

0 F t 


( 47 ) 


(48) 


(49) 


(50) 


(51) 


in which F J is an orthogonal trigonometric function satisfying the condition of zero stress at the 
endpoints, i.e., 

F i (s = s 0 ) = F i (s = si) = 0 (52) 

These functions are formulated as 


Fi(s ) = sin 


iVr(6>(^)-6>(^ 0 )) 

(0(ai)-0(s o )) 


with i = l,7 


(53) 


in which 0 (£q) and ^(^j ) represent the beginning and end angles of the contact region. These angles are 
measured with respect to the global X-axis, and any point in the contact region is identified by the angle 

0(s). 
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The expression for the unknown traction vector at the contact region can be rewritten as 




( 54 ) 


where 


F = [F 0 F, F 2 ■■■ F/] and ={A^ 0 , Aj H) A ( /°} (55) 

From Eqs. (49) and (54), the boundary integral of the product that appears in the 

expression for the total potential energy is obtained as 


J i m c (k ° dr = J A (U) Y T T^\jM a W dr- l A m F T G (kC) A^ dT 
+ \ A (k()T F T t>* m dr 


(56) 


f(*0 


or 


J i^ T c m df - A ,M > C«« a«> - A (t,), 'g<*«A < '> + A (( ° r g<‘' ) (57) 


r(kO 


where 

Q(kO = J pT’ pW 

g{*°= J F r G (W) </r 

f m 

g ( o° = J F r s^r 

p(^0 

Also, the potential energy of the traction vector, , acting on is expressed as 

| t < ‘ )T u<‘>rfr= {u«>V>rfr= Ja<‘» 7 'u<‘>t< t Vr=a« |r f 

p(^) pW p(^) 

where 


(58) 


(59) 


(60) 


( 61 ) 
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Substituting from Eqs. (B3 1), (41), (57), and (61) for the appropriate terms in Eq. (28), the expression for 
the total potential energy is obtained in matrix form as 




k= 1 k= 1 

«<*>-£ 


^=1 


*= I 
K L 


k= 1 


*=1 £=1 
K 

-IT< 

*=1 


a ( *> 2A ( ^ )r g! W) A w +XXA (H)7 'g^ ) 
*=1 £=1 


j: l 


(63) 


*=U=1 


(*) 


in which 


h«>=± f s«V‘>u<‘> T rfr 

(64) 

2 \ 

p(^) 


•h<‘> = J • N «) r n<‘>U«» I 'rfr 

(65) 

p(^) 



where and U^are given by Eqs. (B31) and (B32). The total potential energy in Eq. (63) is 
compacted to its final form 


n = ^a r Ha -*ha + ^-b A + A r Ca - A r f + ACa - A r g, A + A r g 0 -fa 
in which the vectors and matrices are defined by 



( 66 ) 

(67a) 

(67b) 

(67c) 

(67d) 
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J T ={f (1)7 \ t {2)T f (K)T 


(67e) 


v=<|v ,)r , v 2 > r ,..„ v k > t 


(670 


f r ={fd) r , f (2) T fW 7 


(67g) 


zT *(2>t a ( ^> r lwith & ik)T a (k2)T a {kL)T 

So — 180 ’ Sq ’—>So ? with gQ i So ’Sq ’—’Sq 


(67h) 


and 


H 


H 


( 1 ) 


H 


( 2 ) 


H 


( K ) 


(68a) 


,0) 


»(2) 


JL) 


(68b) 


C = 


c (l) 


c (2) 


■<K) 


(68c) 


C = 


’(1) 


’( 2 ) 


’(K) 


withC (t) = 


£(*I) T £(k2) T ... £(*£/ 


(68d) 


and 
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(68e) 





with = 





The minimization of the total potential energy, i.e., Sn = 0 , leads to the following equilibrium equations 


H 0 C T C T r a j ff + *h 

0 b 0 g[ Ml 0 

C 0 0 0 | A [ I f 

C g| 0 0 J l A J l go 


(69) 


Solving Eq. (69) for the unknown vectors permits the calculation of the stress and displacement 
components in each laminate, bolt deflections, and forces exerted by the bolts. Because the angles d A 
and 6 b defining the contact region are unknown, these equations are solved for their assumed values 
until convergence is achieved through an iterative scheme. This procedure is explained in Appendix D. 


4. NUMERICAL RESULTS 

Three different types of load transfer through bolted joints are considered in order to validate the present 
analysis. The first configuration is a pin-loaded square plate considered by Ireman et al. (1993). The 
second configuration is a single-lap joint with four bolts investigated by Xiong and Poon (1994). The 
third configuration is a joint of dissimilar bars with a single bolt under different uniform temperature 
changes analyzed previously by Gatewood (1957). Then, the capability of the present analysis is 
demonstrated by considering two different double-lap joint configurations: one with three bolts under 
mechanical loading only, uniform temperature change only, and their combination, and the other with 
seven bolts under mechanical loading only. 

A pin-loaded square plate was investigated through detailed finite element analysis by Ireman et al. 
(1993). As shown in Fig. 5, the geometry of the plate is defined by the parameters W = 24 mm, D - 6 
mm, 5= 0.021 mm, and thickness t = 3.046 mm. Three different sets of laminate material properties 
considered in the analysis are given in Table 1. The bolt was assumed to be rigid and the applied load, P, 
was taken as 5483 N. In the present analysis, the number of terms, N, retained in the series representation 
of the complex potential functions, <J>^^(z£*^) and z , is taken as 80. The number of terms 

retained in the series representation of the reaction tractions, X^and denoted by J and 7, 
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Table 1. Material properties of laminates. 


Case 

E x , GPa 

E y , Gpa 

Gyy, GPa 

v xv 

A 

99.2 

35.5 

8.5 

0.24 

B 

35.5 

99.2 

8.5 

0.08 

C 

51.5 

51.5 

19.3 

0.33 


*X 


Fig. 5. Pin-loaded plate configuration. 


respectively, is taken as 24. Figures 6-8 shows the favorable comparison of the normalized radial and 
tangential stresses obtained through the finite element analysis reported by Ireman et al. (1993) with those 
of the present calculations. The radial and tangential stress components are normalized with respect to the 
applied bearing stress, -P/(d t ) , and the applied stress, P/(W t ) , respectively. 

The geometry and loading conditions of a single-lap joint of aluminum and composite plates with 
four bolts are shown in Fig. 9. As shown in this figure, the geometry of the joint is defined by the 
dimensions of h ^ =0.31 in, h ^ =0.117 in, D = 0.3125 in, W= 3.125 in, 5 = 1.25 in, e = 0.9375 in, i - 
2.75 in. There exists no bolt-hole clearance, 8 - 0 . The composite lay-up is [(45707-4570 o )2/0790 o ]$ 
with lamina properties of E L = 18.5xl0 6 psi , E T =1.9xl0 6 psi , Gu =0.85 xlO 6 psi , and v L j =0.3. 
In obtaining the present predictions, the series representations of the functions are truncated at 
N = 30 and 7=7=5. Although the entire geometry of the lap joint is considered in the present 
analysis, only the predictions for the bolt loads exerted by Bolt 1 and Bolt 3 are presented due to the 
presence of symmetry. As compared in Table 2, the present analysis predictions are in remarkable 
agreement with those calculated by Xiong and Poon (1994). 

The geometry of the steel and aluminum plates connected with a single bolt is shown in Fig. 10. Each 
plate is subjected to a different uniform temperature change. The bolt is assumed to be rigid, with no 
clearance between the bolt and the bolt holes. The material properties and the temperature change for 
each bar is presented in Table 3. In the present analysis, the bars are assumed as narrow plates whose 
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Normolized circumferential 
stress, a 00 /(P/(Wt)) 



-90 -60 -30 0 30 60 90 


Angle, 0 



Angle, 0 


Fig. 7. Comparison of normalized stress around hole for pin-loaded plate for laminate type B 



Normolized circumferential 
stress, a w /(P/(Wt)) 



Angle, 0 



Angle, 9 


Fig. 8. Comparison of normalized stress around hole for pin-loaded plate for laminate type C. 













Table 3. Material properties and temperature changes 



Material 

E, psi 

V 

aj° F 

> 

o 

•n 

Plate 1 

Steel 

3x 10 7 

0.3 

6.5 x 10 6 

500 

Plate 2 

Aluminum 

10 7 

0.3 

12 x 10' 6 

100 


geometry is specified by W = 0.5 in, e - 3.5 in, s = 0.5, D = 0.2 in, and h ^ =h ^ =0.06 in. The present 
analysis predictions are obtained by considering N = 30 and J = I = 5 in the series representations of the 
functions. The present analysis predicts the force exerted by the bolt to be 442 lb. The bolt force 
obtained by Gatewood (1957) according to the strength of materials approach is 461 lb. The difference of 
4.1% in bolt load prediction is possibly due to the fact that the strength of materials approach does not 
account for the presence of the bolt hole. 

The geometry of a double-lap joint with three bolts is shown in Fig. 11. The geometry of the plates 
all having the same dimensions are described by W = 50 mm, s = 25 mm, e = 12.5 mm, D = 6 mm, h = 5 
mm, 4 mm, and 5 = 0.005 mm. The material properties for the steel and aluminum plates, 

respectively, are: E s = 200 GPa, v s = 0.3, a s = 1 1.7 x 10' 6 /°F, E A = 70 GPa, v A = 0.3, and a A = 23.0 x 10‘ 6 
/°F. All of the bolts in this joint are assumed to be rigid. For the first case of mechanical loading only, 
the applied tensile stress, cr 9 is 200 N/mm. In the second loading case, the plates are under uniform 
temperature change only. Accordingly, the steel plates are free from temperature change, A T s = 0°C, 
while the aluminum plate is subject to A T A = 125°C. The third case represents a combination of the 
mechanical and thermal loads, with <7 - 200 AVmm, AT S = 0°C, and A 7^ = 125°C. The number of terms 
retained in the series representation of the functions in the present analysis is specified by 
N = 30 and J = / = 5 . 

Due to the presence of symmetry, only the results concerning Bolt 2 and Bolt 3 are presented in Figs. 
12-14 and Table 4. The tangential and radial stresses around the bolt holes, shown in Figs. 12-14, 
correspond to mechanical loading, uniform temperature change, and their combination, respectively. As 
observed in these figures, the location and extent of the contact region, as well as the magnitudes of the 
bearing stresses, change significantly. The bolt load distributions for each of these loading cases are 
given in Table 4. As expected, the force exerted by Bolt 2 alters its direction when the nature of the 
loading changes from mechanical to thermal. This behavior is caused by the thermal expansion of the 
aluminum plate constrained with three bolts. Consequently, under thermal loading, the bolt forces with 
their directions toward the middle of the plate resist the thermal expansion of the plate. 
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Fig. 13. Stresses around bolt holes in alloy plate due to thermal loading only. 



Fig. 14. Stresses around bolt holes in alloy plate due to combined thermal-mechanical loading. 
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Table 4. Bolt load distribution in a double-lap joint of steel plates and an aluminum 


plate with three bolts. 




Bolt load, N 

Mechanical Loading 

Thermal Loading 

Combined Loading 

Bolt # 2 

F x 

-2313.99 

6906.94 

4383.20 

F y 

-54.74 

6597.65 

6372.89 

Bolt # 3 

F x 

-5374.18 

-13815.46 

-18769.19 

F^_ 

0.00 

0.00 

0.00 


Under combined mechanical and thermal loading, the forces exerted by Bolt 2 and Bolt 3 are 
approximately equal to the sum of the forces acting on these bolts in the cases of mechanical and thermal 
loads. This shows the effect of the non-linearity arising from the contact analysis. 

A complex double-lap joint with seven bolts, shown in Fig. 15, is subjected to a tensile loading. The 
outer plates are steel and the inner plate is a composite. The bolt-hole clearance is specified as 1% of the 
hole diameter. The composite plate is made of graphite and fiberglass with material properties 
E x = 4.7xl0 5 6 psi, =4.75xl0 6 psi , Gu = 1.2xl0 6 psi , and =0.24. Steel properties for the 
plates and bolts are taken as E = 30.10 6 psi and v = 0.3 . The applied load is 70,000 lbs, corresponding 
to a - 10,000 psi. 

In the calculation of the results, the number of terms retained in the series representation of the 
functions in the present analysis is specified by TV = 30 and 7 = 7=7. Because of the presence of 
symmetry, only the stress distributions around Bolts 1, 2, 4, and 5 are shown in Fig. 16. The forces 
exerted by these bolts are given in Table 5. As expected, the bolts in the first row share more of the load 
than the bolts in the second row. Also, the outer bolts share more of the load than the inner bolts. 
Accordingly, Bolt 1 and Bolt 3 are the most highly loaded bolts for this configuration. 

5. CONCLUSIONS 

In this study, a new approach based on a complex potential theory in conjunction with a variational 
formulation has been introduced for the thermo-elastic contact analysis of a general bolted-joint 
configuration containing multiple laminates joined by multiple bolts. The total potential energy of the 
joint is formulated by using a solution in the form of a complex potential series that automatically 
satisfies the stress equilibrium equations and compatibility conditions, thus avoiding the necessity to 
perform area integrations and resulting in boundary integral expressions for the strain energy of the 
laminates. The total potential energy also includes the strain energy of the bolts based on a shear 
deformable beam theory. 
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Table 5. Bolt load distribution in a double-lap joint of 
steel plates and a composite laminate with 
seven bolts. 


Bolt 

F x , lb 

Fv, lb 

Fy IF' applied 

1 

-1181.5 

11983.7 

0.1712 

2 

0.0 

11340.9 

0.1620 

l ' 4 

2.5 

9947.8 

0.1421 

| 5 

238.3 

7427.6 

0.1061 


In order to capture high gradient variations of stresses near the free or bolted holes, the stress field is 
defined as the superposition of complex potential series originating from each hole. Hence, not only are 
continuous stress and displacement fields obtained but the modeling of the entire joint is also simplified 
considerably. By only entering boundary information, hole size and locations, and the number of terms to 
be used in complex and other series, the solution provides all the stress, displacement, and contact force 
distributions at any point in the joint. 

Contact between the bolts and the laminates is established by enforcing displacement continuity along 
the contact region between the bolts and the plates. This is established by incorporating the work done by 
the unknown contact forces over the contact displacements into the total potential energy expression. The 
contact displacements are defined by constraint equations that take into account the gap between the bolts 
and plates. The contact forces are assumed in the form of trigonometric series that satisfy stress-free 
conditions at the ends of the contact regions. 

Since the contact regions are unknown a priori , an iterative scheme is adopted in order to determine 
the beginning and end angles of the contact regions. Starting with an initial guess, the system matrix is 
generated to solve for unknown plate and bolt displacements and contact forces simultaneously. The 
simultaneous determination of bolt displacements and contact forces, along with the plate displacements, 
is a unique feature of the present formulation. A new guess is then obtained by monitoring the stress 
distribution along the contact regions. The iterative scheme is continued until a configuration for contact 
regions is reached where all the contact forces become compressive. 

The validation problems show excellent agreement of the present formulation against those reported 
by other investigators. The pin-loaded panel in the first validation problem provides a comparison of 
contact angles and contact force distribution. For all laminate configurations, remarkable agreement is 
obtained between the present analysis and the refined finite element solution. In the case of a single-lap 
configuration subjected to thermal loading, the strength of material solution is available. As expected, the 
present analysis achieves the right bolt load and compares well with the strength of material solution if 
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full contact is assumed around the bolt. In the case of a single- and a double-lap joint containing multiple 
holes, the present analysis captures the correct load distribution shared by each bolt. 

The versatility of the present formulation has been demonstrated by solving a double-lap joint 
configuration containing three bolts and two laminates with different material properties, and the joint is 
subjected to thermal, mechanical, and thermo-mechanical loadings. All of these cases are solved 
assuming variable contact regions. Therefore, the rule of superposition is invalid for these problems since 
the contact regions are changing. This can be clearly observed in the tabulated results where the 
summation of the bolt load distribution arising from thermal and mechanical loads is significantly 
different from those arising from the combined loading. 
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APPENDIX A 

The bolt stiffness matrix is derived based on the Timoshenko’s zeroth-order shear deformable beam 
theory. The cross-section of a bolt connecting K laminates (regions) is shown in Fig. Al. The bolt 
number is denoted by i and the regions are numbered from bottom to top in sequential order. The bolt 
has a circular uniform cross-section, A ^ moment of inertia, /^, and Young’s and shear moduli, E \ and 
Gi , respectively. The nuts at the ends of the bolt are assumed to represent clamped boundary conditions, 
thus preventing rotations but creating reaction moments at the ends of the bolt. 

The bolt is subjected to forces arising from the contact between the bolt and the laminates. Because 
of the variation in laminate thickness and stiffness, these forces exerted by the laminates vary through the 
length of the bolt. Because of the variable contact forces along the length of the bolt, the large in-plane 
bolt stiffness compared to those of the laminates and the small ratio of the bolt diameter to its length, the 
most suitable and accurate representation of the bolt can be achieved by discretizing the bolt into small 
Timoshenko beam elements connected at nodal points, as shown in Fig. Al. 





L 


K I 


L, 


L, 


Fig. Al. A close view of the section in the vicinity of a bolt, and the discretization 
of the bolt into Timoshenko beam elements. 


Final Report NAG-1 -2052 


Page 32 





The bolt discretization is based on the most effective locations of the contact forces. Thus, two of the 
nodes are selected at the top and bottom ends of the bolt in order to obtain the largest in-plane bolt 
deflection. The contact forces exerted on the bolt by the top and bottom laminates are assigned to these 
end nodal points. The remaining nodal points are chosen at the intersections of mid-planes of the inner 
laminates and the bolt longitudinal axis, as shown in Fig. Al. Hence, the contact forces exerted on the 
bolt by the inner laminates are assigned to these intermediate nodal points. The nodal deflections and 
rotations along the length of the bolt permit the determination of the bolt deflection at any point along the 
bolt by utilizing interpolation functions. 

Because the bolt material is homogeneous and isotropic and its cross-section is circular, the moment 
of inertia of the bolt is homogeneous on the (x ,y) plane. Hence, there exists no coupling between the 
deformations of the bolt on the (*,z) and (y,z) planes, leading to the uncoupled stiffness matrices of the 
bolt associated with the (x,z) and (y,z) planes. However, their forms are identical since the stiffness 
properties of the bolt are homogeneous on the(x,y) plane. Therefore, the derivation of the stiffness 
matrix associated with th e(jt,z) plane applies to the stiffness matrix associated with the (y,z) plane. 

As shown in Fig. Al, the bolt connecting K laminates is modeled by (K-l) number of beam 
elements and K nodes. Each node is assigned a deflection, A^=A^(z,-), and a rotation, 
= (p^\zi) , with the subscript i = l,K representing the node numbers. The positive directions of the 
deflections and rotations are shown in Fig. Al. Also, the length of each beam element is denoted by 
with / = 1, K . Based on the geometry and material properties of the bolt, the strain energy arising 
from the bolt deformation associated with the (x, z) plane can be written as 

= I>f } (AD 

i=l 

The strain energy of each element is expressed as 



u u 


in which the strains K^P and are based on Timoshenko’s zeroth-order shear deformation theory, 
and they are defined as 


K 


m _ . 


2aO‘Z) 


d A 


zz 


and y%P = 


d A 


m 


dz z dz 

In Eq. (A2), A ^ denotes the corrected area of the bolt and is defined as 


m 

X 


(A3) 
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Aft -c A c 


(A4) 


in which c 2 represents the shear correction factor, which is a correction to the strain energy due to 
uniform transverse shear deformations. 

The displacement and rotation field of the bolt is represented by piecewise continuous interpolation 
functions. These functions are defined individually over each element as 


A^\z) = H l (s)A\Y+H 2 (s)A^ M) +H 3 (s)^-f- + H 4 (s) 


AO 


(O 


d A 


(0 


d a ( 5> 


XI 

dz 


x(M) 


dz 


^(z) = N 1 (s)^Un 2 (s)/^ +1) + N 3 (s)42 


(A5a) 


(A5b) 


where the subscript m denotes the mid-point of the element. The variable s = z-Zj is a local coordinate 
system for the i * element. The interpolation functions Hj(s) (j = 1,.., 4) represent the cubic Hermitian 
polynomials defined as 




3 s 2 2s 3 

Z| + L? 


H 2 (s) = 


3 s 2 2 s 3 


L 2 


L i 

„ . , 2s 1 s 3 

H 3 (s) = s — — +- j 
Li I 2 

h 4 (s)=~+ 4 

Li I 2 

Also, the functions Nj(s ) (j = 1,..,3) represent the Lagrangian shape functions defined as 

A/,(s) = (^-Z,/2)(i-Z 1 )/[f|/2] 

N 3 (s) = s(L i -s)/[ll/ 4] 


(A6a) 

(A6b) 

(A6c) 

(A6d) 


(A7a) 

(A7b) 

(A7c) 


In order to express and 0' x u> in terms of the nodal unknowns defined at the end points of the 
element, two successive steps are performed. In the first step, constraint of uniform shearing strain along 
the beam element is enforced into the kinematic field, thus leading to 


*<«‘0 
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dyW_ = d^j d<p { p 

dz dz 1 dz 


(A8) 


Substituting from Eqs. (A3) and (A5) into Eq. (A8) and grouping the terms as coefficients of z° and z , 
this constraint equation produces two algebraic equations of the form 


4 + 2 


,(0 


Lj dz 1^ dz 


(0 


.A a ^+A a (o +±/v + ±m) 

; 2 " + Q *(<+0 + ^ i. 9 x(i+\) Vxm 


(A9a) 


6 dtfg 6 d^' +l = _ 12 (o 12 (o 

Z 2 dz Lj dz Z| ' " 




3 xi 


zf z? 


“*0’+D ' 72^' ' 72 ^(.+D 


8 <*<*> 


(A9b) 


These two equations are then solved for dA xi ; /dz and dA^. + ^/<fe as 


- 1 a(Q | 1 a(Q + 5 .(/)_i /0 _lM) 

dz Xl Z,- *(< + D 6*” 6^ (i+1) 3*”" 


(A 10a) 


d A 


(0 

*2+2.*— L A ( ^+— -!*<*>+ _!*<*) 

dz Xl Jt(,+1) 6™ 6^ (,+1) y xm 


(A 10b) 


Substituting from Eq. (A 10) into Eq. (A5a) and rearranging the terms, the in-plane deflection 
component, A^ , is obtained as 


A<« = 


1-i 




U 


aW+IaW + 

a xi + . a *(i+l) + 


r 5s 3s 2 2s 3 ^ 

+ 


6 2 314 




f 2 ~ 

S S 25 
+ 


3 7 


6 2 1^ 314 


2 

i J 


Cl) + 


^ 2s 2s 2 4s 3 V, 

— + — <b\ 

6 ° ' 


(All) 


V 


^ 3Z?J 


(0 

xm 


Thus, the total number of unknowns in the expressions for the deflection and rotation is reduced from 7 to 
5. Substituting from Eqs. (All) and (A5) for A^ and ^^in the strain definitions, Eq. (A3), and 
carrying out the integrations in Eq. (A2), analytically result in the following strain energy expression for 
the z th element of the ^ th bolt in matrix form 


r (io] 

T 

Ki 


1 J i( ) 


l q x2 \ 



D 11 

b (*> 

°12 

Jd) T 

h m 

t>i2 

"22 


ic 

i«2\ 


(A 12) 


where 
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A ( o A (o jo jo } 

A ' ; A t(t+1) ™ fr(i+l)j 


J l O-JO 

1x2 ~ Txm 


(A 13a) 
(A 13b) 
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(A13d) 


L (/0 _ 16£^// , AG e A f( Li 
°22 ~ + " 


3A 


(A13e) 


A further reduction of the total number of unknowns from 5 to 4 is achieved in the next step. The 
slope defined at the mid-point of the beam, <p^ , can be condensed out of the strain energy expression, 
Eq. (A 12), by employing the static condensation procedure. In the absence of nodal forces at the mid- 
point of the element, the first variation of the strain energy with respect to q^ yields the equilibrium 
equation at the mid-point as 

b lf r q?f ) +b 22 0x2 =° ^ I 4 ) 


Solving for (i.e., <p^ ) in the above equation and substituting into Eq. (A 12), the strain energy 
expression reduces to 


W? 0 4# )r b ( ' 0 q <«> 


in which q^ is identical to , and is defined as 


(A 15) 
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(A 17a) 


(A 17b) 


(A 17c) 


Substituting from Eq. (15) into Eq. (Al), the strain energy of the i * bolt is expressed as 


;= l 1 


(A 18) 


As mentioned previously, the presence of nuts prevents the bolt from rotating at the end points; thus, it 
resembles clamped-type boundary conditions, requiring that =0 and Invoking this 

condition into Eq. (A18), the strain energy expression is modified as 


U 


(0 


= 2 q < 


m T b 'do 


'(10 


, K-2 
Z 1=2 


<lx 


MO 




1)0 


h '((*- DO n '((AT-I)0 


(A19) 


where the vectors and matrices with a prime are defined as 
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(A20d) 


The terms on the right-hand side of Eq. (A19) can be rearranged such that all the unknown deflection and 
rotation components and the coefficients of the element stiffness matrices are assembled in large vectors 
and matrices. Hence, the strain energy expression is compacted to 
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°24 +0 13 


b (0 - 

u A A — 


.(20 

5 23 


b (f) = 


6 ( tf > +fc <20 
^44 + °33 


b (2C) 

°34 


u(20 

i4 


b m +b m 

°24 ±a \3 


b oo 

04 


h ((K- 2)0 
°23 


b ((K- 2)1) y 
0 24 *°\3 


(( K-m 


h ((K- 1)0 
°23 


b (2{) 

°34 


b (20 A 30 
^44 + ^33 


, ( 30 

^34 


" 7 “ 

I 

l - 

- 4 - 


((K-3)C) 

? 34 


,((K-3)0 > h ((K -2)C) 
°44 + "33 


A(K- 2)0 
D 34 


A(K-2)0 

°34 


(A22d) 




((AT- 2 )r) +/b ((^-!)0 


44 


33 


The bolt is subjected to contact forces that are assigned to the nodal points and the reaction moments due 
to clamped support at the end points of the bolt. 

None of the internal nodal points, where the nodal rotations are active, is subjected to external 
moments. Therefore, it is appropriate to reduce the total number of unknowns by statically condensing 
out the rotational components of the bolt. Because no external moments are acting on the internal nodal 
points, the first variation of the strain energy with respect to the vector q^ 1 must vanish, thereby yielding 
the moment equilibrium equations in matrix form as 


b 


S ^A )+b #^=° 


(A23) 


Solving 


in the above equation and substitution into Eq. (A21), after rearranging the terms, lead to 


</<'> =J-Ai f > r b< f >A<P 


where the vector is identical to q' V and the matrix b^ ' is defined as 


( 0 . 




w 1 kW 7 K «’r 1 K (o 


b [° =b\ 1 ’ --b 


AA 2 # 


(A24) 


(A25) 


As mentioned previously, the stiffness properties of the bolt in the (x,z) and (y,z) planes are identical 
because the bolt material is isotropic and its cross-section is circular. 
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Thus, the strain energy of the bolt in the (y,z) plane can be expressed as 


{/ (f) = I A (0 7 ' b (f) A (0 
u y 2 y y y 


in which the vector contains 


A ( 0 ' = / A (0 a (0 ... A ( 0 l 

*>' I vl J2 ^yK\ 


and the matrix is identical to , i.e., 


(0 


h (0 _ b <*> -Ib w V )_1 b W 

D y ~ d aa 2 D # 


The total strain energy of the bolt becomes 




or, in a more compact form, 


(/<«=! 


K>] 

T 

i 

O 

/--s 

w H 

pO 

1 


k°l 

k°j 


. 0 


K’j 


(A26) 


(A27a) 


(A28b) 


(A29) 


(A30) 


The analytical derivation of Eq. (A30) depends on the number of plates connected with the ^ th bolt. For 
example, the analytic derivation of matrix (=b^ ) for bolts used in single- and double-lap joints is 
obtained as 


with 



(A31) 


* 12/ij 

~ LjG(Af£ + l2E(I(Li 


(A32) 


for a single-lap joint, and 
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b'° =■><;> = 


t oo _ 

"ll 


b (u) 

_41 


■ (If) , ,(20 
"44 + "22 


Sym. 


.,(10 ,.,(10 , ,,(2fk 
.(If) "14 ("34 + "12 ) 

13 .(10 ..,<20 

"44 +"22 



fc (10,(2f) 

"14 "23 
, (If) , (20 
"44__+"22__ 

,JlO~ .(2fk,(2f) 

i (If) . "34 + "12 '"23 

"l3 +' 


1,(20 _ 

"33 


h (U) +b {2C) 

"44 ___22 

n (20 2 


23 


^44 + °22 


(A33) 


with ^ and & ; y 2 ^ defined by Eq. (A16) for a double-lap joint. If a bolt connects more than 3 
laminates, the analytic derivation would be too lengthy to present herein. For this case, numerical 
calculation of the matrices becomes more appropriate. 


APPENDIX B 

The strain energy given as 


"“’4 J N $ e $ M - J 

Aik) ,<*) 


(Bl) 


can be rewritten in the form 


{/<*>=! J jN<« r V«^ 

2 Ak) A (k) 


(B2) 


in which 


N^ r ={^,<,<} 


.(*)' _ r J) Jk) Jk ). 
& v °jcjc ’ c >^y » C jcy J 


* (*)' _rM*> V*> 

fc \ C Xt , Cyy , / 


Substituting for the stress resultant vector, , given as 

N (*) =A (*) e W 

leads to 


(B3) 
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(B4) 


U (k) =- J e ik)T A. (k V k) dA-- Je^A « V»<M or 

2 a (*) 2 >«(*) 


Its first variation becomes 


su m = f {a'V'-a' 1 ' v*>} 7 '& ( *>« 

d(*> 


Utilizing the kinematic relations 


o ( * ) = i/ u ( U .„(*) A 


and invoking Eqs. (3) and (6) into Eq. (B5) result in 


d(*) 


Before applying the divergence theorem, this expression can be rewritten as 

' P *<*) 


*(*) 


, f Is \ 

Noting that N a p ^ — 0 and applying the divergence theorem yield 

«/«>= /{^-W<g]»^i4‘>dr- | NtyfSjPdA 


<k) 


Finally, this expression can be recast as 


5U {k) =- J N ( JJ pSu {k) dA+ J *t {k) Su (k) dT 

A (k) ’ f(k) 

L with a,fi = x,y and p,7] = n 

+ j '$>&$> d r + X j '.fSufdr 

p(^) ^=1 p(£-0 


in which 


with 


,(*) 


and 


#>*#> +4*v«f =*<*> *«?> +/?>&?> =#>&.<» 


(B5) 


(B6) 


(B7) 


(B8) 


(B9) 


,s (BIO) 


(Bll) 
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APPENDIX C 


As introduced by Lekhnitskii (1968), the displacement field, and«^, and the in-plane stress 
resultants, N* k \ Ny k \ and , satisfying the equilibrium equations and compatibility condition for the 
k 0 ' region (laminate) can be written in terms of arbitrary complex potential functions 


and 



(Cl) 


(C2) 


The complex constants, p r * and q rk (r = 1,2) , are given by 


Prk = a \\f4k + a \2 ~ a l 6 W ( C 3) 

c !rk =a \2 > Mrk +a 22^rk ~ a 26 (^ 4 ) 


in which the complex parameters ju^ and fok are the roots of the characteristic equation derived by 
Lekhnitskii (1968), 


(*)„4 M) . „(Ak,,2 (*)„ _ n tr-i z\ 

a l\ t*rk & rk + & a \2 +l3 66 ^rk ^ a 66 ^rk + a 22 (^5) 

where a^\ with (/, j = 1,2,6), are the compliance coefficients of the A region (laminate). 

For the k ^ region containing multiple bolt-holes, the complex potential functions, <J>^ and <p ^ , can 
be expressed by superposing the complex potentials, and <p ^ , for individual bolt holes as 

L (k) 

4 ,(t ) ( 2 « ) ) = ^ 4> (K > ( z (« ) ) (C6) 

f = l 

and 



(C7) 


where is the number of bolt holes in the laminate and and with 

r = l,2, are the complex functions for the k th region containing the £** hole located at (X The 
complex parameters and are defined as 
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Z «> = X +/ i, ( y and 


in which X and Y are the global coordinates and jc^ and are the coordinates associated with the 
bolt hole. The explicit form of these complex potential functions can be expressed in the form 

t (#*r (cs) 

n=-N ' 

n*0 

and 



(C9) 


in which are complex unknown coefficients and map the f 1 * 1 circular hole in the k * 

laminate to a unit circle in the mapped plane, thus permitting Laurent series representations to 
approximate the field variables. The prime denotes differentiation with respect The mapping 

functions introduced by Bowie (1956) are in the form 


em = 



a kt (1 " ‘Mrk ) 


(CIO) 


in which a k/ > is the radius of the £ ,h hole in the k' h region and i = -J - I . The sign of the square-root term 
is chosen so that |^*^|> 1 . Inverting the mapping function provides 

= <C11) 

in which 

= SK {l - i/irt ) , s «fl = ^.( 1+1>rt ) (C12) 

The displacement components, and can be rewritten in terms of real vector quantities as 


i k2, 

f=\ r — 1 M V 


,U 


m 


f=l r =ln=-N 


x(rn) '^y(rn) 


AM) 


(C13) 
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(Cl 4) 


L {k) 2 N 




' s m T K «o 7 t-Ko 7 l.om 

M=1 r=l/z=-/A 


J XT(r/i) ’^>y(rn) ,3 ry(rn) l“™ 


where the real vectors are given by 

T 


<2) ={ 2Re 




U ( / ( ^) 7 ={2Re[ 9r , (W) <„],-2Im[^ (W) <,]] 


; (W) ' 
’ xx(rn) 


= {2Re 


,,2 (k£)* 1 o T „J,,2 (Jfcf),** 

/V* VrrtJ*“ 2Im [/^rA: Vr/i 


]) 


S »(m) r ={ 2Re [- <if) ^ ].-2ta[- <W V; ]} 


o(*0 j 

°*y(r/«) 


= {2Re 


Mrk (W Vr«] -2Im[// r * ( ^V™]] 


,(*0' 


= { 2 Re [ ( (i U) cc rn ] , -2 Im[ {k£) a m ] } 


The functions <J>^* and are defined as 


and 


®!T(£‘ 0 H£‘ 0 )" 


with 


<Pi 


(k£)* _ 


dz 


m 


n(£ ke) \ 

,U(^U )Z1 L 

X m J ,.MO 


O), 


MO -MO . -V 

— ’r “ r 


m 


(£“>) 


Terms arising from the expansion of , S an< * a r//' for r ranging from 1 to . 

contained in the following vectors: 


it(^0 ^ _ J|j(^0 ^ tt(^0 ^ 
u cr(n) _ 1 ^a(ln) ’ u a'(2n) 


(Cl 5a) 
(Cl 5b) 
(Cl 6a) 
(Cl 6b) 
(Cl 6c) 
(Cl 7) 

(Cl 8) 
(Cl 9) 

(C20) 
can be 

(C21) 
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o(M) ' 
°afi(n) 


= S 


: (k() 


’#/?(!«) ’^a/?(2n) 


(kC) 


(C22) 




- L (kt) 1 

' a In 


,a 


(*0 7 

2« 


(C23) 


The terms arising from the expansion of these vectors for n ranging from -N to N are contained in the 
vectors defined as 


II 

D 

/ijCAO T |j(*0 T 

| ^ai-N) ’' J a(-N+ 1) ’■ 

T 

" ,VJ a{N) 

s {k i )T 

o m T MkC) T 

°afi(-N) '*aP{-N+\) »■ 
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■■'°aP(N) 

a tkC)T 

_L(«) 7 .(«) r 

| (-N) ’“(-AM-l) 

( ^) r l 

■’V) J 


(C24) 

(C25) 

(C26) 


with (a, P = x, y) . Thus, using Eqs. (C24)-(C25), the series expansions for displacements and resultant 
stresses in Eqs. (C13) and (C14), respectively, can be rewritten as 


,(k) 


,<*)- 


'a 




Akt) 


a v 


£=l 


(C27) 


L (k) 

™ a /3 Zu a 

^=1 

with (a = x, y) . These expressions can be recast as 


jik) 


u ( *) = £ U ( ^ a w 
f=l 


(C28) 


(C29) 


by defining the following vectors 


L (0 

s (k) = Y l S (k ° T a m 
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(*) J 




N 1 


( k)‘ 


= { 


N ( Jt\N (k) 




(C30) 
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u m =[i4 W) uf 


i(kt) _ 


c(&0 c(kO 

^ XT °>'>’ °X)' 


Finally, these equations can be expressed in a more compact form as 


u<*>=U ( *> a W 


where 


N (*) =s (*) r a (*> 


u ( *> r = 


u 


(* 1 ) 


r v tn> T ... uH *’) 7 


Aky 


s (kl) T s (k2) T ___ g 




Ak) T _ 


Jk\) T Jkl) T 

a , a v , a 




(C31) 

(C32) 


APPENDIX D 

The iterative scheme for solving the system of algebraic equations given in Eq. (69) begins with the initial 
estimates of 9^ and 9^ , shown in Fig. Dl, defining the contact region. These angles are measured 
in the counterclockwise direction from the x ^ axis of the local coordinate system, (;c^,y^). The 
initial estimate in most cases does not represent the true contact region, for which the radial stresses are 
all compressive a fP <0 on F ( ^ , and the fact that the start and the end points of the contact region 
have zero radial stresses, cr^ ) (9 A ) = 0 and crj.P(9 s ) = 0. 

As shown in Fig. Dl, three distinct cases exist, depending on whether the radial stress near the 
start angle, 9^ , is larger or smaller than the true value of the start angle or equal to the true value, 9 A , 
of the contact region, f ' ^ . The initial guess of the starting angle, 9 ^ , is smaller than its true value, 
9 a , if the compressive radial stresses change sign and become tensile near the start point of the contact 
region, 9^ . The initial guess of the starting angle, 9^ , is larger than its true value, 9 A , if the 
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Fig. Dl. The behavior of radial stress near the start point of a contact region: 
(a) the start angle is too small; (b) the start angle is too large; 
(c) the start angle is correct. 


compressive radial stresses do not change sign and remain compressive near the start point of the contact 
region. The initial estimate of the end angle, 6 ^ , is larger than its true value, 0 B , if the compressive 
stresses change sign and become tensile near the end angle. If the compressive stresses do not change 
sign, then the initial guess, 6 ^ , is smaller than its true value. 

During the iteration process, an initially guessed contact region defined by 6^ and 6 ^ converges 
to the true contact region defined by 0 A and 0 B through incrementally changing the values of the start 
and end angles. The increment is forced to decrease each time the true value of the start or end point is 
passed, and the direction of angle change is altered. The convergence of the iterative process is achieved 
when the incremental value of the angle reaches a pre-defined value. 

In order to avoid the case of radial stresses having zero values at the start or end points of the contact 
region but with tensile and compressive stresses along the contact region, two auxiliary points are 
considered as shown in Fig. D2. These points are located inside the contact region near the start and end 
points of the contact region. In order to achieve convergence for the contact region, the radial stresses at 
these two auxiliary points must be compressive. 
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